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LONG-TERM  GOALS 

The  long  term  objective  of  our  research  for  the  “High  Resolution  Air-Sea  Interaction”  (HRES) 
Departmental  Research  Initiative  (DRI)  is  to  identify  the  couplings  between  large  wave  events, 
winds,  and  currents  in  the  surface  layer  of  the  marine  boundary  layers.  Turbulence  resolving 
large  eddy  simulations  (LESs)  and  direct  numerical  simulations  (DNSs)  of  the  marine  atmo¬ 
spheric  boundary  layer  (MABL)  in  the  presence  of  time  and  space  varying  wave  fields  will  be 
the  main  tools  used  to  elucidate  wind-wave-current  interactions.  A  suite  of  turbulence  simula¬ 
tions  over  realistic  seas  using  idealized  and  observed  pressure  gradients  will  be  carried  out  to 
compliment  the  field  observations  collected  in  moderate  to  high  winds.  The  database  of  simu¬ 
lations  will  be  used  to  generate  statistical  moments,  interrogated  for  coherent  structures,  and 
ultimately  used  to  compare  with  HRES  observations. 


OBJECTIVES 

Our  near  term  goals  are:  1)  participate  in  the  planning  process  for  the  HRES  field  campaign; 
and  2)  construct  an  LES  code  applicable  to  the  HRES  high  wind  regime.  In  order  to  accomplish 
the  latter  goal  we  are  improving  the  parallelization  of  our  base  LES  code  and  developing  an  al¬ 
gorithm  to  allow  simulations  of  turbulent  winds  over  nearly  arbitrary  3-D  wave  fields. 


APPROACH 

We  plan  on  investigating  interactions  between  the  MABL  and  the  connecting  air-sea  interface 
using  both  LES  and  DNS.  The  waves  will  be  externally  imposed:  (1)  based  on  well  established 
empirical  wave  spectra;  or  (2)  ultimately  provided  by  direct  observations  of  the  sea  surface  from 
field  campaigns.  The  main  technical  advance  is  the  development  of  a  computational  tool  that 
allows  for  nearly  arbitrary  3-D  wave  fields,  i.e.,  the  sea  surface  elevation  h  =  h(x,y,t)  as  a  sur¬ 
face  boundary  condition.  The  computational  method  will  allow  time  and  space  varying  surface 
conditions  over  a  range  of  wave  scales  (9(10)m  or  larger. 
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WORK  COMPLETED 


Meetings:  We  attended  a  PI  meeting  in  Monterey  CA  that  focused  on  refining  elements  of  the 
HRES  field  campaign.  The  discussion  concentrated  on  the  atmospheric  measurements  to  be  col¬ 
lected  from  aircraft  and  from  R/V  FLIP  with  particular  attention  paid  to  the  surface  layer  pres¬ 
sure  measurements.  A  summary  of  the  discussion  is  available  from  C.  Friehe  (U.C.  Irvine).  We 
also  presented  a  summary  of  our  modeling  results  from  the  HRES  and  ITOP  (Impacts  of  Ty¬ 
phoons  on  the  Western  Pacific  Ocean)  departmental  research  initiatives  at  the  ONR  program 
review  (June,  2009)  in  Chicago. 

Publication:  Last  year  we  wrote  an  article  for  Annual  Review  of  Fluid  Mechanics  (see  Sullivan 
&  McWilliams,  2010).  Our  contribution  discusses  the  coupling  processes  between  surface  gravity 
waves  and  adjacent  winds  and  currents  in  the  marine  boundary  layers.  Wind-wave  and  wave- 
current  interactions  are  important  as  these  processes  modulate  the  exchanges  of  momentum, 
heat,  and  gases  between  the  atmosphere  and  ocean.  Atmospheric  processes  we  discuss  are  wind- 
driven  waves,  wave-driven  winds,  steep  waves  and  drag  laws.  In  the  ocean  boundary  layer,  we 
summarize  the  currently  accepted  views  of  wave-current  interactions,  the  modeling  of  Langmuir 
circulations,  impacts  of  breaking  waves,  and  the  implications  of  wave-current  interactions  for 
mixing  at  high  winds. 

Algorithmic  Developments:  During  the  past  year  we  began  adapting  our  parallel  LES  code 
(Sullivan  &  Patton,  2008)  to  accommodate  a  general  3D  time  varying  wavy  surface.  This  was 
accomplished  in  the  following  stages:  a)  Develop  the  governing  equations  for  an  atmospheric 
PBL  in  curvilinear  time-dependent  surface-following  coordinates;  b)  Build  and  test  a  module 
capable  of  generating  a  3D  wavy  surface;  and  c)  Implement  and  test  a)  and  b)  in  our  “flat  bot¬ 
tom”  parallel  LES  code  for  the  atmospheric  PBL. 

Introducing  a  3D  moving  bottom  boundary  adds  significant  complexity  compared  to  our  previ¬ 
ous  boundary-layer  simulations  with  a  single  monochromatic  wave  (Sullivan  et  al.  2007).  Here, 
we  briefly  outline  a  few  key  elements  of  the  problem  formulation  and  computational  implemen¬ 
tation  for  the  general  wavy  bottom  problem  with  emphasis  on  the  treatment  of  advection  and 
pressure  [further  details  are  available  in  Sullivan  (2010)].  In  the  present  formulation  the  ad¬ 
vection  of  momentum  and  scalars,  and  also  subgrid-scale  energy  in  the  case  of  LES,  (ptq,  9,  e) 
is  written  in  flux  conservative  form  using  a  contravariant  velocity  that  takes  into  account  grid 
movement.  A  generic  transport  equation  for  variable  if  in  time-dependent  surface  following  co¬ 
ordinates  takes  the  form 


I© +  4[^-^)1  =  7  (1) 

where  £  =  &  =  (£,  ?/,  C)  are  computational  coordinates,  Uj  is  the  contravariant  flux  velocity,  zt 
is  the  grid  speed  (i.e.,  the  speed  at  which  vertical  grid  lines  translate),  J  is  the  Jacobian  con¬ 
necting  physical  and  computational  spaces  x  =>■  £,  and  the  right  hand  side  1Z  includes  the  appro¬ 
priate  combination  of  pressure,  buoyancy,  viscous,  Coriolis,  and  diffusion  terms.  Since  the  grid 
is  constrained  to  translate  vertically  the  grid  motion  only  contributes  vertical  flux  ifzt.  Also, 
included  in  our  equation  set  is  the  geometric  or  space  conservation  law  (SCL)  (Thomas  &  Lom¬ 
bard,  1979)  which  relates  the  change  of  an  elementary  control  volume  to  the  co-ordinate  frame 
velocity.  Inspection  of  the  advective  term  in  (l)1  shows  that  with  constant  if  and  uniform  veloc¬ 
ity  the  SCL  is 

1  The  space  conservation  law  follows  most  naturally  from  the  continuity  equation  written  for  a  variable  den¬ 
sity  fluid. 
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With  an  externally  imposed  surface  wavefield,  either  J  or  zt  can  be  specified  at  future  time 
steps.  Then  the  SCL  is  satisfied  by  computing  the  counterpart  zt  or  J  according  to  the  partic¬ 
ular  discretization  of  (2).  In  our  implementation,  we  choose  to  compute  zt  based  on  knowledge 
of  J  at  forward  time  steps. 

Our  boundary-layer  model,  using  DNS  or  LES,  adopts  the  incompressible  Boussinesq  assump¬ 
tion  and  thus  an  elliptic  Poisson  equation  for  a  pressure  variable  p  must  be  solved  in  order  to 
satisfy  the  continuity  equation.  This  is  the  major  computational  issue  in  modeling  incompress¬ 
ible  flows  above  terrain.  The  pressure  Poisson  equation  for  a  surface  following  grid,  developed 
from  the  continuity  equation  expressed  in  terms  of  contravariant  flux  velocities,  dUi/d^i  =  0,  is 
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where  S  is  the  divergence  of  the  right  hand  side  of  the  momentum  equations  and  (/ x ,  / y ,  (z) 
are  metrics  arising  from  the  coordinate  transformation.  Direct  solution  of  this  elliptic  equation 
is  not  readily  possible  using  the  standard  method  of  2D  Fourier  transformation  in  horizontal 
planes  followed  by  a  tridiagonal  matrix  solve  for  each  pair  of  horizontal  wavenumbers  ( kx ,  ky ) 

( e.g .,  Sullivan  et  al.  2000).  This  is  due  to  the  variable  coefficients  and  mixed  derivatives  result¬ 
ing  from  the  grid  transformation.  Instead,  we  enforce  mass  conservation  by  iteratively  solving 


(n+l)  _ 


=  V^p(n)  + 


A  t'y 


V  •  U(n) , 


(4) 


for  the  pressure  p.  In  (4),  n  is  the  iteration  index,  At  is  the  time  step  and  7  is  a  weight  associ¬ 
ated  with  the  time  integration.  W2Dp  is  a  diagonal  pre-conditioning  matrix  that  can  be  inverted 
using  standard  methods.  Tests  indicate  that  this  iterative  scheme  is  robust  and  converges  the 
solution  for  p  to  machinery  accuracy  with  roughly  n  =  30  iterations  for  waves  of  modest  slope 
ak  ~  0.1.  More  demanding  tests  with  steep  waves  (and  hills)  are  planned. 

A  wavy-bottom  DNS/LES  boundary- layer  model  is  completed  by  specifying  the  wavefield  h  = 
h(x,  y ,  t )  and  appropriate  boundary  conditions.  As  a  first  step,  we  build  h(x,  y,  t )  from  a  linear 
combination  of  simple  plane  waves.  The  amplitudes  and  phases  of  the  components  are  picked  to 
mimic  properties  of  the  wave  field,  e.g.,  short  or  long  crested  waves,  wave  propagation  direction, 
or  measured  wave  height  spectra  {e.g.,  Pierson  and  Moskowitz,  1964).  This  flexibility  allows  us 
to  conduct  idealized  process  studies  with  a  known  wavy  surface  but  also  provides  the  framework 
to  ingest  measured  wavefields  from  the  HRES  field  campaign.  Periodic  boundary  conditions  are 
used  in  the  horizontal  (/,  p)  directions  as  is  customary  for  boundary  layer  simulations.  At  the 
air-water  interface,  the  atmospheric  motions  are  set  equal  to  the  the  water  orbital  velocities.  At 
the  present  time,  the  wave  motions  are  assumed  to  be  irrotational. 

The  numerical  method  used  to  solve  the  transport  equations  for  momentum,  scalars  and  subgrid- 
scale  energy,  the  pressure  Poisson  equation,  and  the  space  conservation  law  are  broadly  similar 
to  those  described  in  Sullivan  et  al.  (2000,  2007).  The  spatial  differencing  is  pseudospectral  along 
computational  coordinate  lines  (/,  p)  and  second-order  finite  difference  along  the  vertical  coordi¬ 
nate  line  (.  A  third-order  Runge-Kutta  time  stepping  scheme  operating  with  a  fixed  CFL  num¬ 
ber  is  used  to  advance  the  equations  in  time.  We  emphasize  that  identical  numerics  are  used  to 
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solve  the  momentum  and  space  conservation  laws  in  order  to  avoid  false  generation  of  sources 
and  sinks  due  to  the  grid  movement  (e.g.,  Ferziger  and  Peric,  2002).  A  fractional  step  projection 
method  is  used  to  enforce  mass  conservation  on  the  contravariant  flux  velocities.  The  most  ex¬ 
pensive  step  of  the  solution  procedure  is  the  iterative  solution  of  the  pressure  Poisson  equation 

(4). 

The  parallel  implementation  of  the  wavy-bottom  algorithm  outlined  above  utilizes  a  2D  domain 
decomposition  similar  to  our  “flat  bottom”  turbulence  simulation  code.  Each  processor  operates 
on  three-dimensional  “bricks”  sub-sampled  in  x,  y,  or  2  directions.  Brick-to-brick  communica¬ 
tion  is  a  combination  of  transposes  and  ghost  point  exchange.  To  preserve  pseudospectral  differ¬ 
encing  in  the  horizontal  directions  a  custom  parallel  matrix  transpose  was  designed  and  imple¬ 
mented.  Finally,  a  new  pressure  Poisson  solver  was  coded  to  match  the  domain  decomposition 
[further  details  are  given  in  Sullivan  &  Patton,  (2008),  Sullivan  (2010)]. 


RESULTS 

The  computational  algorithm  and  component  parts  described  in  the  WORK  COMPLETED 
section  are  currently  being  implemented  and  evaluated  for  idealized  cases.  Figures  1  and  2  are 
typical  snapshots  of  instantaneous  wavefields  generated  by  our  virtual  wave  algorithm.  In  each 
figure,  the  wave  amplitudes  are  chosen  to  match  the  Pierson-Moskowitz  wind-wave  equilibrium 
spectrum  for  winds  of  15  ms”1.  In  figure  1  the  waves  propagate  with  no  directional  dependence 
while  in  figure  2  long  crested  waves  propagate  across  the  grid.  In  these  examples,  the  space  and 
time  scales  of  the  wavefields  lie  in  the  resolved  scale  regime  of  an  atmospheric  LES.  Figures  3 
and  4  are  illustrative  examples  of  LES  of  atmospheric  boundary  layers  with  complex  surface 
layers.  In  figure  (3),  we  show  the  near  surface  boundary-layer  wind  field  around  an  isolated  3D 
hill  of  height  50  m.  The  winds  are  driven  by  a  mixture  of  surface  heating  Q*  =  0.05  K  ms”1 
and  geostrophic  winds  ( Ug,Vg )  =  (5,0)  ms”1.  The  boundary-layer  depth  is  (9(1000)  m.  Notice 
the  speedup  of  the  winds  near  the  crest  and  sides  of  the  hill.  In  figure  (4)  we  show  contours  of 
the  horizontal  wind  at  a  height  of  6  m  above  a  two-component  propagating  surface  wave  field. 
These  trial  simulations  demonstrate  that  the  new  LES  code  is  functional  and  capable  of  han¬ 
dling  complex  lower  boundaries.  In  the  future  we  plan  a  detailed  validation  of  the  wavy  bottom 
simulation  code  using  the  rough-hill  wind  tunnel  measurements  described  by  Gong  et  al.  (1996) 
and  the  wind- wave  tank  observations  described  by  Veron  et  al.  (2007).  These  are  demanding 
test  cases  as  they  include  regions  of  extensive  flow  separation. 


IMPACT  /APPLICATIONS 


The  computational  tools  developed  and  the  database  of  numerical  solutions  generated  will  aide 
in  the  interpretation  of  the  observations  gathered  during  the  HRES  field  campaigns.  In  addition 
idealized  process  studies  performed  with  the  simulations  have  the  potential  to  improve  parame- 
terizations  of  surface  drag  under  high  wind  conditions  in  large  scale  models. 


TRANSITIONS  &  RELATED  PROJECTS 

We  are  currently  engaged  in  analyzing  data  collected  during  the  Ocean  Horizontal  Array  Tur¬ 
bulence  Study  (OHATS)  and  the  Coupled  Boundary  Layers  Air-Sea  Transfer  (CBLAST)  field 
campaigns.  These  are  joint  efforts  between  NCAR,  and  numerous  university  investigators.  Also 
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the  present  work  has  links  to  the  ONR  DRI  on  the  impact  of  typhoons  in  the  Western  Pacific 
Ocean  (1TOP). 
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Figure  1:  An  example  of  a  virtual  wave  field  generated  by  linear  combinations  of  plane  waves. 
The  amplitudes  of  the  individual  components  are  picked  to  match  a  Pierson-Moskowitz  wind- 
wave  equilibrium  spectrum  at  a  wind  speed  of  15  m  s_1.  The  wavefield  is  sampled  at  a  reso¬ 
lution  of  Ax  =  Ay  =  2.9  m  and  propagates  from  left  to  right.  The  wave  height  amplitude 

h  =  h(x,  y,  t )  shown  in  the  color  bar  is  in  units  of  meters. 


Figure  2:  An  example  of  a  wave  field  propagating  across  the  x  —  y  computational  grid  with  a 
directional  distribution  function  chosen  to  emphasize  long  crested  waves.  The  wind  speed  and 
grid  resolution  are  identical  to  those  in  figure  1. 
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Figure  3:  The  turbulent  wind  field  (u—  component)  around  and  above  an  isolated  hill  in  a 
mixed  convective-shear  driven  atmospheric  boundary  layer  -  surface  heat  flux  Q*  =  0.05  K  m 

s^1  and  geostrophic  wind  Ug  —  5  m  s-1.  The  3D  cosine  shaped  hill  is  50  m  tall.  Notice  the  wind 
speedup  over  the  crest  and  side  of  the  hill  with  u/Ug  >  1.2.  The  color  bar  at  the  top  of  the 

figure  is  in  units  of  m  s^1.  This  simulation  employs  256  x  256  x  128  gridpoints  and  demonstrates 
that  the  LES  algorithm  is  capable  of  simulating  turbulent  flow  above  general  topography. 


Figure  4:  The  horizontal  u— component  of  the  wind  field  at  a  height  of  6  m  above  a  moving 
two  component  wave  field  produced  by  large-eddy  simulation.  Pockets  of  high  and  low  winds 
are  concentrated  around  the  crests  and  troughs  of  the  wavefield.  The  color  bar  at  the  top  of  the 

figure  is  in  units  of  m  s-1.  The  simulation  employs  256  x  256  x  128  gridpoints. 

8 


